#include <iostream>
#include <vector>
#include <math.h>

#include "HEALPIX_Data.h"
#include "fitsio.h"

using namespace std;

int compare_pixels (const void *a, const void *b) {

  healpix_t *pix1;
  healpix_t *pix2;


  long in1;
  long in2;
  long temp;

  pix1 = (healpix_t *) a;
  pix2 = (healpix_t *) b;

  in1 = pix1->pixel;
  in2 = pix2->pixel;
  temp = in1 - in2;
  if (temp > 0)
    return 1;
  else if (temp < 0)
    return -1;
  else
    return 0;
}
int HEALPIX_Data::Get_Sky_Val(float sky_ra, float sky_dec, vector<float> &values) {

	/* first get the healpix theta and phi */
	long hpix;
	long ns = nside;
	double phi = sky_ra;
	double theta = M_PI/2.0 - sky_dec;
	void *result = NULL;	
	healpix_t key;
	ang2pix_ring(ns,theta,phi,&hpix);
	
	key.pixel = hpix;
	result = bsearch((void *) &key, (void *) &data[0], data.size(), sizeof(healpix_t), compare_pixels);

	if (result == NULL) {
		return EXIT_FAILURE;
	}
	else {
		healpix_t *found = (healpix_t *) result;
		values[0] = found->value[0];
		values[1] = found->value[1];
		values[2] = found->value[2];
		values[3] = found->value[3];
		return EXIT_SUCCESS;
	}
}	
	 
int HEALPIX_Data::Get_Pix_Val(int pix, vector<float> &values) {

	healpix_t *the_pixel;
	values.resize(4);
	if (verbose)
		cout << "HEALPIX_Data::Get_Pix_Val Entering the Get_Pix_Val method in the derived calss HEALPIX_Data" << endl;
	
	the_pixel = &this->data[pix];

	values[0] = the_pixel->value[0];
	values[1] = the_pixel->value[1];
	values[2] = the_pixel->value[2];
	values[3] = the_pixel->value[3];
	
	return EXIT_SUCCESS;

}


int HEALPIX_Data::Get_Pix_Vertices(int pix, vector<float> &point, vector<float> &vertices) {

/* -----------------------------------------------------------------------------
 *
 *  This member function is based on a slight modification of earlier work by others.
 *
 *  The slight modification is the extension of this routine to the calculation of
 *  pixel vertices by Stephen Ord (2010). And the transliteration to C++ and the 
 *  transformation from a standalone function to a class member function. 
 *  This is Released as part of the MWAVIS package under 
 *  the GNU GPL v2.
 *
 *  This is a slight modification of a file distributed as part of
 *  HEALPix which is	
 *
 *  Copyright (C) 1997-2008 Krzysztof M. Gorski, Eric Hivon, 
 *                          Benjamin D. Wandelt, Anthony J. Banday, 
 *                          Matthias Bartelmann, 
 *                          Reza Ansari & Kenneth M. Ganga 
 *
 *
 *  HEALPix is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  HEALPix is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with HEALPix; if not, write to the Free Software
 *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 *
 *  For more information about HEALPix see http://healpix.jpl.nasa.gov
 *
 *----------------------------------------------------------------------------- */

  int nl2, nl4, npix, ncap, iring, iphi, ip, ipix1;
  double  fact1, fact2, fodd, hip, fihip;
  double PI=M_PI;
  double HALFPI = M_PI/2.;
  double z,phi;
  // vertices specific
  double phi_nv, phi_wv, phi_sv, phi_ev;
  double z_nv, z_sv;
  double hdelta_phi;
  double theta;
  int iphi_mod, iphi_rat,ipix;

  char err_string[1024];

  //      PARAMETER (pi     = 3.1415926535897932384626434d0)
  //      parameter (ns_max = 8192) ! 2^13 : largest nside available
  
  int ns_max=8192;
  // locus
  healpix_t *the_pixel;
  if (verbose)
	  cout << "HEALPIX_Data::Get_Pix_Vertices--In the Get_Pix_Vertices method in the derived calss HEALPIX_Data" << endl;
 
	
  the_pixel = &this->data[pix];
  ipix = the_pixel->pixel;
	

  try { 
	  if( nside<1 || nside>ns_max ) {
		  sprintf(err_string, "%s (%d): nside out of range: %d\n", __FILE__, __LINE__, nside);
		  throw(err_string);
	  }
	  npix = 12*nside*nside;      // ! total number of points
	  if( ipix<0 || ipix>npix-1 ) {
		  sprintf(err_string, "%s (%d): ipix out of range: %d\n", __FILE__, __LINE__, ipix);
		  throw(err_string);	
	  }
  
	  ipix1 = ipix + 1; // in {1, npix}
	  nl2 = 2*nside;
	  nl4 = 4*nside;
	  ncap = 2*nside*(nside-1);// ! points in each polar cap, =0 for nside =1
	  fact1 = 1.5*nside;
	  fact2 = 3.0*nside*nside;
  
	  if( ipix1 <= ncap ) {  //! North Polar cap -------------
    
		  hip   = ipix1/2.;
		  fihip = floor(hip);
		  iring = (int)floor( sqrt( hip - sqrt(fihip) ) ) + 1;// ! counted from North pole
		  iphi  = ipix1 - 2*iring*(iring - 1);
    
		  z = 1. - iring*iring / fact2 ;

		  theta = acos(1. - (pow(iring,2)/(3. * (pow(nside,2)))));
		  phi   = (1.*iphi - 0.5) * PI/(2.*iring);

		  phi_nv = phi;
		  phi_sv = phi;


		  // Vertices calculation
		  hdelta_phi = PI/(4.0*iring);   //  half pixel width
		  z_nv = 1.0 - pow((iring-1),2.) / fact2;
		  z_sv = 1.0 - pow((iring+1),2.) / fact2;
		  iphi_mod = (iphi-1) % iring;  // ! in {0,1,... iring-1}
		  iphi_rat = (iphi-1) / iring; //     ! in {0,1,2,3}
		  if (iring > 1) 
			  phi_nv = HALFPI * (iphi_rat +  iphi_mod   /double(iring-1));
		  phi_sv = HALFPI * (iphi_rat + (iphi_mod+1)/double(iring+1));
	

	  }
	  else if( ipix1 <= nl2*(5*nside+1) ) {//then ! Equatorial region ------
    
		  ip    = ipix1 - ncap - 1;
		  iring = (int)floor( ip / nl4 ) + nside;// ! counted from North pole
		  iphi  = (int)fmod(ip,nl4) + 1;
    
		  fodd  = 0.5 * (1 + fmod((double)(iring+nside),2));//  ! 1 if iring+nside is odd, 1/2 otherwise
		  z = (nl2 - iring) / fact1;
		  theta = acos( (nl2-iring)/(1.5 * nside));
	
		  phi   = (1.*iphi - fodd) * M_PI /(2.*nside);
		
		  // vertices calculation
		  hdelta_phi = PI/(4.0*nside); //  ! half pixel width
		  phi_nv = phi;
		  phi_sv = phi;
		  z_nv = (nl2 - iring +1) / fact1;
		  z_sv = (nl2 - iring -1) / fact1;
		  if (iring == nside) { //! northern transition
			  z_nv = 1.0 - pow((nside-1),2.0) / fact2;
			  iphi_mod = (iphi-1) % nside; // ! in {0,1,... nside-1}
			  iphi_rat = (iphi-1) / nside; //     ! in {0,1,2,3}
			  if (nside > 1) 
				  phi_nv = HALFPI * (iphi_rat +  iphi_mod   /double(nside-1));
		  }
		  else if (iring == 3*nside) { // ! southern transition
			  z_sv = -1.0 + ((nside-1)*(nside-1)) / fact2;
			  iphi_mod = (iphi-1) % nside; //! in {0,1,... iring-1}
			  iphi_rat = (iphi-1) / nside; //     ! in {0,1,2,3}
			  if (nside > 1) 
				  phi_sv = HALFPI * (iphi_rat +  iphi_mod   /double(nside-1));
		  }

	  }
	  else {//! South Polar cap -----------------------------------
    
		  ip    = npix - ipix1 +1;
		  hip   = ip/2.;
		  /* bug corrige floor instead of 1.* */
		  fihip = floor(hip);
		  iring = (int)floor( sqrt( hip - sqrt(fihip) )) + 1;//     ! counted from South pole 
		  iphi  = (int) (4.*iring  + 1 - (ip  - 2.*iring*(iring-1)));
		  theta = acos(-1. + (iring)*(iring)/(3.*nside*nside));  // !is this correct??
		  z = -1. + iring*iring / fact2 ;
		  phi   = (1.*iphi+0.5) * M_PI/(2.*(iring));
		  // Vertices calculation
		  phi_nv = phi;
		  phi_sv = phi;
		  hdelta_phi = PI/(4.0*iring); //  ! half pixel width
		  z_nv = -1.0 + pow((iring+1),2) / fact2;
		  z_sv = -1.0 + pow((iring-1),2) / fact2;

		  iphi_mod =(iphi-1) % iring ; //! in {0,1,... iring-1}
		  iphi_rat = (iphi-1) / iring ; //    ! in {0,1,2,3}

		  phi_nv = HALFPI * (iphi_rat + (iphi_mod+1)/double(iring+1));
		  if (iring > 1) phi_sv = HALFPI * (iphi_rat +  iphi_mod   /double(iring-1));


	  }


	/* 

	we have all the info now to build everything - we will assume clockwise pixel ordering
	N,E,S,W - we will present them in theta, phi.

	There is a slight correction as healpix uses latitude and we need co-latitude (M_PI/2. - theta)

	*/
	
	  point[0] = phi;
	  point[1] = theta;

	// north vertex

	vertices[0] = phi;
	vertices[1] = acos(z_nv);

	// east vertex
	
	phi_ev      = phi + hdelta_phi;
	vertices[2] = phi_ev;
	vertices[3] = theta;

	// south vertex

	vertices[4] = phi;
	vertices[5] = acos(z_sv);

	// west vertex

	phi_wv = phi - hdelta_phi;
	vertices[6] = phi_wv;
	vertices[7] = theta;
	
	return EXIT_SUCCESS;


  }
  catch (const char *e) {
	cerr << e << endl;
	return EXIT_FAILURE;
  }
}
		
int HEALPIX_Data::load() {

	cout << "HEALPIX_Data::load Entering the load method in the derived calss HEALPIX_Data" << endl;
	return EXIT_SUCCESS;
}


HEALPIX_Data::HEALPIX_Data (char *filename) {

	/*! Now this has to load everything in -
	We know the format a FITS BINARY TABLE - so we need to use some fitsio calls
	to get at the data. This comes from the RTS routine load_mwahpx_from_file

	*/

	float nullval=0;;
	int anynul=0;
	fitsfile *fptr=NULL;
	int status=0,nfound=0,hdutype=0;
	long ns=0,naxis=0;
	long *naxes=NULL;
	long nrow=0,nrows=0,pcount=0;
	int tfields=0;
   	char extname[64];
	
	
	cout << "HEALPIX_Data::HEALPIX_Data(char *filename) Overloaded Constructor" << endl;
	nullval = HEALPIX_NULLVAL;

	try {
 
		if (fits_open_file(&fptr,filename,0,&status)) {
			throw("fits_open_file has failed\n");
		}

		/* currently the binary table is in the next hdu */
   
		if (fits_movrel_hdu(fptr,1,&hdutype,&status)) throw("Cannot move to next HDU\n");
   
		if (fits_get_hdu_type(fptr, &hdutype, &status)) throw("Cannot get HDUTYPE\n");
   
		if (!(hdutype == BINARY_TBL)) {
			throw("ERROR: file does not contain MWAHPX data. Or wrong HDU tested. Current version of MWAHPX expects BINARY TABLE in the second HDU\n");
		}
		else {
    
			if (fits_read_key_lng(fptr, "NAXIS", &naxis, 0x0, &status)) {
				throw("fits_read_key has failed on NAXIS\n");
			}
    
			naxes = (long *) malloc(naxis*sizeof(long));

			if (fits_read_keys_lng(fptr, "NAXIS", 1, (int)naxis, naxes, &nfound,&status)) {
				throw("fits_read_key has failed on NAXIS\n");
			}

			if (nfound > 1)
				nrow = naxes[1];

			if (!nrow) {
				throw("ERROR: Could not determine image size\n");
			}

			/* Number of pixels per side of each base-resolution pixel. */

			if (fits_read_key_lng(fptr, "NSIDE", &ns, 0x0, &status)) {
				throw("NSIDE not found\n");
			}
			cout << "I think NSIDE is " << ns << endl;

			if (fits_read_btblhdr(fptr, 32, &nrows, &tfields, NULL, NULL, NULL, extname, &pcount, &status)) {
   				throw("Cannot read binary table\n"); 
			}
			fprintf(stdout,"EXTENSION %s\n",extname);
			if (strncmp(extname, "WEIGHTED", 8) == 0) {
				fprintf(stdout,"STATUS: Weighted!\n");
				throw("Cannot yet read a weighted file\n");

			}

			/*! Need to allocate some space for the HEALPIX Data */

			this->data.resize(nrow);
			this->nside = ns;
		

			for (int irow=0; irow < nrow ; irow++) {
				if (verbose)
					cout << "Reading Row: " << irow << "\r" ;

				if (fits_read_col(fptr, TINT, 1, irow+1, 1,1,
							&nullval, &data[irow].pixel, &anynul, &status)) {
					throw("Failed column read\n");
				}

				if (fits_read_col(fptr,TFLOAT,2,irow+1,1,1,
							&nullval, &data[irow].value[0], &anynul, &status)) {
					throw("Failed column read\n");
				}

				if (fits_read_col(fptr,TFLOAT,3,irow+1,1,1,
							&nullval, &data[irow].value[1], &anynul, &status)) {
					throw("Failed column read\n");
				}


				if (fits_read_col(fptr,TFLOAT,4,irow+1,1,1,
							&nullval, &data[irow].value[2], &anynul, &status)) {
					throw("Failed column read\n");
				}

				if (fits_read_col(fptr,TFLOAT,5,irow+1,1,1,
							&nullval, &data[irow].value[3], &anynul, &status)) {
					throw("Failed column read\n");
				}

				if (fits_read_col(fptr,TFLOAT,6,irow+1,1,1,
							&nullval, &data[irow].weight[0], &anynul, &status)) {
					throw("Failed column read\n");
				}
			
	
			}

			cout << "\nRead " << nrow << " rows in total" << endl;
			ndata=nrows;
		}
	}
	catch(char const *e) {
		fits_report_error(stderr,status);
		cerr << e << endl;
	}

}

void HEALPIX_Data::Get_Val_Limits(vector<float> &limits) {

	healpix_t *current = NULL;

	limits.resize(8);

	vector<float> I;
	vector<float> Q;
	vector<float> U;
	vector<float> V;

	I.resize(ndata);
	Q.resize(ndata);
	U.resize(ndata);
	V.resize(ndata);

	for (int pix=0;pix<ndata;pix++) {

		current = &this->data[pix];
		
		I[pix] = current->value[0];
		Q[pix] = current->value[1];
		U[pix] = current->value[2];
		V[pix] = current->value[3];
		
	}

	stable_sort(I.begin(),I.end());
	limits[0]=*I.begin();
	limits[1]=*(I.end() -1);
	stable_sort(Q.begin(),Q.end());
	limits[2]=*Q.begin();
	limits[3]=*(Q.end()-1);

	stable_sort(U.begin(),U.end());
	limits[4]=*U.begin();
	limits[5]=*(U.end()-1);

	stable_sort(V.begin(),V.end());
	limits[6]=*V.begin();
	limits[7]=*(V.end()-1);

	I.resize(0);
	Q.resize(0);
	U.resize(0);
	V.resize(0);

}
	
HEALPIX_Data::HEALPIX_Data () {

	cout << "HEALPIX_Data::HEALPIX_Data() Default Constructor" << endl;
}
HEALPIX_Data::~HEALPIX_Data () {

	cout << "HEALPIX_Data::~HEALPIX_Data() Default Destructor" << endl;
}
	
